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BISTABILITY,  BASINS  OF  ATTRACTION,  AND  PREDICTABILITY 
IN  A  FORCED  MASS-REACTION  MODEL 


1.  INTRODUCTION 

/  Bistable  phenomena  can  occur  in  many  physical,  chemical,  and  biological  models  of  natural 
phenomena  [e.g.,  1,  2,  3,  4,  5,  6,^7]?  An  important  subset  of  problems  exhibiting  bistability  consists  of 
those  models  employing  mass-reaction  kinetics:  Crucial  to  understanding  any  mass-reaction  model  is  a 
knowledge  of  the  parameter  corresponding  to  the  rate  of  contact  between  two  or  more  species.  In  cer¬ 
tain  applications,  the  contact  rate  may  be  time  dependent,  and  in  fact,  periodic.  For  example,  the  pro¬ 
cess  of  temporally  increasing  and  decreasing  the  solar  intensity  respectively  changes  the  probability  of 
contact  between  two  reacting  species  in  the  atmosphere'  [8,  9,  10],  In  addition  to  perturbing  reactants 
in  the  atmosphere,  periodic  forcing  of  contact  rates  plays  an  important  role  in  modelling  recurrent  epi¬ 
demic  outbreaksTlfT7)-  It  is  important  to  note  that  both  physical  and  biological  phenomena  exhibit 

r - 1 

oscillations  which  are  longer  than  the  forcing  period,  or  not  periodic  at  all  [10,  121.  Furthermore,  it  is 
not  uncommon  for  periodically  forced  differential  equation  models  to  exhibit  two  or  more  stable 
subharmonic  solutions  for  a  given  set  of  parameters  [e.g.,  2,  7,  131.  The  question  we  consider  here  is, 
how  well  can  one  predict  the  asymptotic  final  state  given  initial  conditions  having  finite  precision  for  a 
problem  that  exhibits  two  different  stable  periodic  orbits. 

■fivs 

In  particularise  consider^  simple  mass-reaction  model  with  a  periodically  forced  contact  rate. 
The  following  is  shown  numerically: 

ytd  j^There  exist  parameter  values  for  which  the  model  exhibits  at  least  two  distinct  stable  subhar¬ 
monic  periodic  orbits. 

J 

h  fiD^jThe  basins  of  attraction  of  each  orbit  can  be  very  complicated,  thus  affecting  final  state 
predictability  as  a  function  of  precision  in  initial  conditions^ 

'Mitfli'prhe  dimension  (capacity)  of  the  basin  boundary  is  estimated  using  the  techniques  in  [14]. 


2.  A  SIMPLE  MASS-REACTION  MODEL 

Let  c,  denote  the  fraction  of  concentrations,  /  -  1,4.  The  model  used  is: 
c{  -  #*(1  -  c()  -  /3 (r)c,c3 

c'i  /3(t)ctCj  —  (a  +  n)Ci  (MR) 

Cj  *  acj  —  (y  +  /i)cj 
ci-  yc3  - 
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where  /jl,  a,  y  are  constants.  Periodicity  in  the  contact  rate  is  incorporated  by  assuming 

0(r)  -  0O(1  +  0i  cos  2 irt),  and  the  forcing  amplitude,  0t,  is  the  parameter  to  be  varied.  Note  that 

* 

JV,  -  1  for  all  time,  and  therefore  it  is  sufficient  to  study  the  first  three  equations  in  (MR),  c4  being 

3 

given  by  1  -  J^c,. 

When  0!  -  0,  equation  (MR)  has  two  steady  states:  (1,  0,  0,  0)  and  -c°  -  (c?,  c$,  c?,  c® ) 
where  c?  -  Q*  +  a)(p.  +  y)/0o a  «■>  l/Q,  c§  -  (41  +  y)ft(Q  —  l)/0oa,  c?  -  n(Q  -  l)/0o-  Parame¬ 
ters  are  chosen  so  that  (1,  0,  0,  0)  is  unstable  and  -c 0  is  stable.  Assume  Q  >  1 ,  fi(Q  —  D  is  small, 
and  introduce  the  change  of  variables:  n  (Q  -  1)  —  «,/*  +  «  “  A^«,  m  +  y  “  A  3/e,  where  e  is  a  small 
parameter  and  A2,  A3  are  positive  constants.  Letting  c,  -  c(°(l  +  x,),  /  —  1,  2,  3  and  using  the  above 
scaling  results  in  the  system: 

x[  -  -€[(•»)  +  0i  cos  2irr)xi  +  0+0!  cos  2irr)x3 

+  0!  cos  2 irt  +  X1X3G  +  0i  cos  2irt)] 

x2  A2(0i  cos  2 irt  +  (xi  +  x3)(l  +  0i  cos  2 irt)  (MS) 

—  x2  +  XiX3(l  +  0i  cos  2irt)/t 
x3  —  A3(x2  —  x3)/e,  where  rj  —  Q/(Q  —  1)>  1. 

The  non-trivial  equilibrium  point  is  moved  to  the  origin,  and  the  eigenvalues  of  the  linearized  vector 
field  at  the  origin  are  given  by  X  ±  (e)  —  tr  ±  iv  +  0(«2),  X3(e)  —  —  (A2  +  A3)/e  +  0(e)  where 
r  -  [A2A3  -  (A2  +  A3)2tj]/2(A2  +  A3)2  <  0,  and  v1  -  A2A3/(A2  +  A3).  Thus  there  is  a  strong  attrac¬ 
tion  onto  a  surface  in  the  direction  corresponding  to  \3,  and  a  weak  attraction  in  which  orbits  slowly 
spiral  into  the  origin.  Since  the  orbits  appear  to  lie  approximately  in  two  dimensions  where 
c2(r)  ~  (4  +  y)c3(t)/a,  we  examine  variations  in  Ci  and  c3.  The  initial  conditions  (ci,  c2,  c3),  are 
restricted  to  satisfying  c2  -  c3(#t  +  y)/a  from  here  on. 

3.  BISTABLE  BEHAVIOR 

For  given  values  of  0i  >  0,  there  exist  periodic  orbits  which  have  period  1  as  well  as  subhar¬ 
monic  periods.  Both  stable  and  unstable  orbits  were  computed  as  a  function  of  0i  using  the  techniques 
in  (15].  The  fixed  values  of  the  parameters  used  in  the  computations  are  listed  in  Table  1.  If  <P(t) 
denotes  the  vector  solution  to  (MS)  of  a  periodic  solution  of  period  n  for  a  given  value  of  0b  define 
the  norm  of  ♦  as 


•  '  *'■>  *  »  *  *  *  »'♦  1**  ,  •  , 
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||4>(f)||2-  (l/n)f*Q(t)  Q(t)dt. 

A  plot  of  ||4>||  as  a  function  of  /3i  is  shown  in  Fig.  1.  Emanating  from  the  stable  steady  state  is  a  small 
amplitude  period  1  orbit.  At  a  critical  value  of  /3(,  the  period  1  orbit  goes  unstable  (dashed  line)  and  a 
stable  period  2  bifurcates  from  the  period  1  branch.  Notice  that  there  also  exist  large  amplitude 
saddle-node  (stable-unstable)  bifurcations  of  periods  3  and  4.  Examples  of  the  period  3  and  period  2 
orbits  are  shown  in  Figs.  2a  and  2b.  From  Fig.  1,  it  is  clear  that  there  exist  at  least  two  stable  periodic 
orbits  for  certain  values  of  Thus  the  model  exhibits  bistability  for  several  ranges  of  forcing  ampli¬ 
tude. 

Table  1  —  Parameter  Values  Used  for  Numerical  Simulation 

#o  -  1575 

H  -  0.02 

a  -  1/0.0279 

X  -  1/0.01 

4.  PREDICTABILITY  AND  DIMENSION  OF  THE  BASINS  OF  ATTRACTION 

When  /31  -  0.1,  Fig.  1  shows  that  there  co-exist  stable  period  1  (SP1)  and  stable  period  3  (SP3) 
orbits.  To  find  the  basins  of  attraction,  we  simplify  the  problem  by  assuming  initial  conditions, 
(cj,  c2,  cj),  are  restricted  by  satisfying  c2  -  (Ox  +  y)/a)e3  as  described  above.  We  pick  a  pair  (cu  c3) 
at  random  using  a  uniform  probability  distribution  and  then  determine  if  the  trajectory  converges  to  an 
SP1  or  SP3  orbit.  If  it  converges  to  SP1,  (cj,  c3)  is  plotted.  Otherwise,  the  trajectory  converges  to  SP3 
and  the  point  is  not  plotted.  (No  other  periodic  orbit  was  observed.)  Figure  3a  illustrated  the  basin  of 
attraction  for  SP1  (dotted)  over  a  given  region.  The  blank  areas  denote  the  basin  for  SP3.  A  uniform 
grid  was  initially  chosen,  and  approximately  24  per  cent  of  the  initial  conditions  converged  to  SP3.  The 
same  procedure  was  performed  for  a  small  region  in  Fig.  3a,  except  that  the  points  were  chosen  at  ran¬ 
dom.  Figure  3b  shows  the  result.  Clearly  the  basins  of  attraction  for  SP1  and  SP3  are  intertwined  in  a 
complex  manner. 


In  [14]  it  is  shown  that  a  discrete  map  in  two  dimensions  possessing  a  complicated  basin  structure 
results  in  an  obstruction  to  predictability.  Here  we  follow  their  procedure  for  measuring  the  predictabil¬ 
ity  of  equation  (MR).  Given  a  random  initial  condition,  (cl5  c3),  we  consider  perturbed  initial  condi¬ 
tions  (cj,  c3  ±  «)  for  a  given  «  >  0.  If  at  least  one  of  the  perturbed  initial  conditions  converges  to  an 
orbit  other  than  the  orbit  passing  through  the  unperturbed  initial  condition,  we  say  (cb  c3)  is  uncertain. 
Four  thousand  random  initial  points  were  chosen,  and  the  fraction,  F,  of  uncertain  initial  points  was 
computed.  The  procedure  was  then  repeated  for  several  values  of  i. 

Figure  4  illustrates  the  scaling  of  F  as  a  function  of  «.  Using  a  linear  least  squares  fit,  it  is  found 
that  F  is  proportional  to  e0'68.  Thus  the  capacity  [16]  of  the  basin  boundary  is  approximately  1.32. 

5.  CONCLUSIONS 

We  have  demonstrated  bistability  for  a  simple  periodically  forced  mass-reaction  model.  Further¬ 
more,  we  have  illustrated  that  the  basins  of  attraction  of  the  periodic  orbits  are  intertwined.  The  effect 
of  the  complicated  basin  structure  is  to  obstruct  the  predictability  of  the  final  state  given  imprecision  in 
initial  conditions.  That  is,  given  the  slope  of  0.68  in  Fig.  4,  to  decrease  the  fraction  of  uncertain  initial 
conditions  in  phase  space  by  an  order  of  magnitude  requires  more  than  an  order  of  magnitude  reduc¬ 
tion  in  uncertainty.  Therefore,  requiring  a  high  degree  of  predictability  can  result  in  the  need  to  per¬ 
form  extraordinary  high  precision  computations. 


Fig.  3b  —  Basins  of  /action  depicting  a  blow-up  of  a  small  region  on  Fig.  3a.  The  initial  conditions 
ere  chosen  at  random.  The  fractal-like  structure  is  evident. 
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Fig.  4  —  Fraction  of  uncertain  points  are  plotted  as  a  function  of  imprecision  in  the  initial  conditions. 
The  dimension  of  the  basin  boundary  was  found  to  be  approximately  1.32. 
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